Levy Flights in External Force Fields: Langevin and Fractional Fokker-Planck 

Equations, and their Solutions 

Sune JespersenB 

Institute of Physics and Astronomy, University of Arhus, 8000 Arhus C, Denmark 

oo ; 

Ralf Metzler! 

School of Chemistry, Tel Aviv University, 69978 Tel Aviv, Israel 
Hans C. Fogedby 

f*S . Institute of Physics and Astronomy, University of Arhus, 8000 Arhus C, Denmark, 

and NORDITA, Blegdamsvej 17, 2100 K0benhavn 0, Denmark 
*0 \ (February 1, 2008) 



> 



00 



We consider Levy flights subject to external force fields. This anomalous transport process is 
described by two approaches, a Langevin equation with Levy noise and the corresponding generalized 
Fokker-Planck equation containing a fractional derivative in space. The cases of free flights, constant 
. force and linear Hookean force are analyzed in detail, and we corroborate our findings with results 

from numerical simulations. We discuss the non-Gibbsian character of the stationary solution for 
the case of the Hookean force, i.e. the deviation from Boltzmann equilibrium for long times. The 
possible connection to Tsallis's g-statistics is studied. 
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l, . I. INTRODUCTION 

S3 

In recent years there has been growing interest in anomalous diffusion in various fields of physics and related 
sciences. In one dimension anomalous diffusion is characterized by a mean square displacement of the form 



((Ax) 2 ) oc 2Df, (1) 



deviating from the linear dependence on time found for Brownian motion 0-|(J. The generalized diffusion constant 
' has the dimension [D] = cm 2 sec -7 . 
7-H , Subdiffusive transport (0 < 7 < 1) is encountered in a diversity of systems, including the charge carrier transport 
in amorphous semiconductors [0J|], N.M.R. diffusometry on percolation structures ||, and the motion of a bead in 
a polymer network |l(J . On fractal structures in general, subdiffusion prevails due to the occurrence of holes of all 
length scales |^]. Examples of enhanced diffusion (7 > 1) include tracer particles in vortex arrays in a rotating flow 
p| , layered velocity fields [jl2j , and Richardson diffusion Jl3| . 

Levy flights are used to model a variety of processes, such as bulk mediated surface diffusion and applications 
in porous glasses and eye lenses ]l5p , transport in micelle systems or heterogeneous rocks [H, special problems in 
reaction dynamics |l7j , Single Molecule Spectroscopy |Q , and even the flight of an albatross |19| . 

Among the different frameworks for describing anomalous diffusion are fractional Brownian motion poj , the contin- 
uous time random walk scheme (JtJ^i], fractional diffusion equations [p2 23 1, generalized Langevin and Fokker-Planck 



equations (FPEs) p4|-|28| , and generalized thermostatistics p9],E30l . Common to all these approaches is the violation of 



the Central Limit Theorem of probability theory |3l|]32|] , and this is achieved either by correlations or by long-tailed 
statistics. Levy statistics J ^.331 is an example of the latter, and has been used extensively to model both enhanced 
and dispersive diffusion |6pl,[33||. The two most fundamental properties of the Levy distributions are the stability 
under addition, following from the Generalized Central Limit Theorem valid for Levy distributions, and the asymp- 
totic power law decay. These features are responsible for the anomalous character of the diffusion processes we have 
in mind. 

A fractional Fokker-Planck equation (FFPE) describing anomalous transport close to thermal equilibrium, was 
presented recently p8| . Since it describes subdiffusion in the force-free case, it involves a strong, i.e. slowly decaying 
memory. In the present paper, we focus on FFPEs which are connected with Levy flights, and are based on the 
following Langevin equation for the coordinate x(t) |24|p5||: 
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Here, m is the mass of the diffusing particle, and 7 denotes a friction coefficient. F(x) is the external force field. For 
simplicity we shall work in one dimension, with obvious modifications in the general case. The noise r)(t) is the source 
of the anomalous behavior. We assume r)(t) to be uncorrelated at different times, and to obey Levy statistics |Q. In 
Fourier space we thus define 



p(k) = J , 



d-qe- ikr >p(ri) = exp(-L>|A;|' i ), 



(3) 



where < < 2. The probability density function (PDF) p(r/) has an asymptotic power-law behavior according 



to p{rj) ~ \r) 



Hj,B2LB5| . Here and in the following, we denote the Fourier transform of a function by using the 



explicit dependence on the wave number k, and analogously u for the Laplace transform of a i-dependent function. 
For the special case fi = 2 in Eq. (§), i.e. for a Gaussian noise, we are led back to the Brownian case. In this case 
2D is the variance of the PDF, but even in the general case the parameter D characterizes the width of the PDF in 
some sense. 

The Langevin equation Eq. (J2]) is a stochastic differential equation. Often it is more convenient to work with the 
deterministic equation for the distribution function, the FPE 
Eq. (ph, we are led to the following FFPE |2lp|]: 



For the power-law noise rj(t) defined through 



F(x)W(x,t) 



7m 



DVW(x,t). 



Here, D denotes the generalized diffusion coefficient with the dimension [D] 
derivative in Eq. fl) is defined through its Fourier-transform pq,B3] 



cm^sec 



(4) 



The Riesz fractional 



V M = - 



d d k 
(2n) d 



\k\' 



(5) 



in d dimensions. Note that in the FFPE Eq. (|J) the first order differential operator acting upon the force term is 
not affected by the introduction of the Levy distribution Eq. (||) , see Ref. [Q where a unifying derivation of FFPEs 
from a generalized master equation is discussed. 

In the following we will consider the FFPE Eq. (Q) for the three cases of the free flight (F — 0), the constant force 
F(x) — Fq, and the Hookean force F(x) — —Ax, comparing to the Brownian case as we go along. We will discuss the 
differences to the sub-diffusive FFPE of Ref. |]28[ , where a fractional operator in time is encountered and the spatial 
part of the standard FPE remains unchanged, as well as the possible connection to Tsallis's g-statistics. Numerical 
simulations corroborate our theoretical findings. We then exemplify the method of solution for the Langevin equation 
Eq. (H) for a linear force with an additional drift term. Before drawing the conclusions, we give some remarks on the 
simulations. Some additional calculations on the nature of the correlation functions are presented in the Appendix. 



II. FREE LEVY FLIGHT 



In this case we have to solve the anomalous diffusion equation 

d 



dt 



W(x,t) = DV^W(x,t). 



(6) 



Fourier transforming Eq. (||) and utilizing the definition of the fractional Riesz operator Eq. (^) we have 



g-W(k,t) - -D\k\»W(k,t) 



with the solution 



W{k,t) 



„-£>t|fe|" 



(7) 



(8) 



demanding the sharp initial condition x(0) = 0, corresponding to W(x,0) = S(x) or W(k, 0) = 1. Comparing to Eq. 
(0) we recognize the characteristic function of the Levy distribution, and we thus find in real space the stable law L M : 



W(x,t) = (Dt)-^L^^ 



77 TT 1 ' 1 



(M/aO, (1,1/2) 
(1,1), (1,1/2) 



(9) 
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In Eq. (||), we have expressed the Levy distribution exactly in terms of Fox' iJ-functions jll],[42]]. This result Eq. 
(^) is expected, due to the stable law nature of the underlying Levy distribution. The asymptotic behavior of the 
propagator W(x, t) can be derived from Eq. (^) and reads 

^0M)~]^T ( 10 ) 

for ^^/[Dt] 3> 1, and thus we encounter a divergence of the mean square displacement at all times: (x 2 {t)) = oo. This 
is intuitively clear due to the occurrence of arbitrarily long jumps in the Levy flight, see Fig. pi Mathematically, the di- 
vergence is evident from Eq. (||) by using the properties of the characteristic function: (x n (t)) — i n d n W(k,t)/dk n \ k=Q . 
In order to extract the scaling form implied by Eq. (0) operationally, one could enclose the "walker" in an imaginary 



growing box (see Sec. VI) 



(x 2 {t)) L ~ / dxx 2 W{x,t) ~ t 2 l» . (11) 

This has been implemented numerically, and as can be seen from Fig. 0, where we have a straight line on a log-log 
plot of (x 2 )l as a function of i, for a fixed /i, the expected power law index 2//i according to Eq. ( pi] ) is found. This 
scaling result {x 2 (t))^ ~ t 2 ^ 1 is not to be confused with the mean square displacement (x 2 ) = J dx x 2 W(x,t) = oo 
for the Levy flight, resulting from Eq. (|l0|). However, for fi> 1, the squared absolute mean 

(\x\) 2 = ( f dx\x\W(x,t)) (12) 



converges, and is proportional to (x 2 (t))i from Eq. (|il|), see the discussion in Ref. f43fl . Fig. || shows this propor- 
tionality for (j, = 1.5 Eij]. In Fig. [|, we graph the slope 2//j as a function of fi for a variety of values for the Levy 
index fi, and we obtain excellent agreement with Eq. ([ll]). In the case /i = 2 we see from Eq. (JsJ) and the results 
of the simulations in Fig. ^ that the usual Brownian behavior is recovered. Especially, we obtain the mean square 
displacement 



„2 



•r 



(t)) = 2Dt. (13) 



The properties of free Levy flights could also be obtained directly from Eq. (||) employing the method of characteristic 
functions. This view also allows us to extract the distribution of speeds, which turns out to be Levy distributed. From 
this likewise follows the mean kinetic energy, and we have for a finite mass m of the walker in the case /i < 2 

(^mv 2 ) = oo. (14) 



III. CONSTANT FORCE: DRIFT AND ACCELERATION 



For a constant force F(x) = F , the FFPE Eq. (0) reads 



dt ' dx \ jm 

Returning to the Fourier domain, we recover the equation 



|Wt) = -#- ( ) +DWW(x,t). (15) 



i-W(k, t) = (-ik^ - D\k\f ) W(k, t), ! I .(i ) 



(17) 



dt ' \ 7TO 
which for the propagator, i.e. W(k,t) with the initial condition W(k, 0) = 1, yields 



W(k,t) = exp -t 



ikl9- + D\k\>* 



This is the same Levy distribution as calculated for the free Levy flight in Eq. d|), but at the translated coordinate 
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W(x,t) = W (x-^,t). (18) 

Here Wq refers to the distribution of the free Levy flight. The displacement of the coordinate is due to the balancing 
of the friction against the imposed constant force, i.e. ^mv = Fq, in the Galilei transformed system x — > x — Fot/^ym]. 
Clearly, the analytical form for the solution in ir-space is still given by Eq. (^]), but now for the translated coordinate. 
Thus Levy flights in a constant force field described by Eq. ( |l5| ) are similar to (anomalous) diffusion in a constant 
velocity field 1 4^j4(| . The reason why we go directly to the steady state described by 7771V = Fo is due to the omission 
of the inertial term in the Langevin equation Eq. (||). By including this term it can be shown that we obtain a 
transient contribution. Thus the diffusion in a force field is only equivalent to diffusion in a constant velocity field 
for large times, according to the discussion in Ref. [55fl, see below. In fact the same situation is encountered even 



for the standard diffusion-advection picture 1 45 4qj4q] 7 For long waiting periods in the random walk picture for the 
subdiffusive model, the external force can only act upon the walker, when it is released from a trap, which leads to the 
sublinear dependence (x(t)) oc t 1 found in Ref. fljsfl , where 7 denotes the power law index of the broad waiting time 
distribution. The Poissonian waiting time distribution, on the other hand, which is typical for the Levy flyer, leads 
to the same behavior Eq. ( ^l|) below [for 1 < fi < 2], as known from the Brownian case, due to the very short waiting 
periods in between jumps. The long-tailed nature of the Levy flight however makes all higher moments infinite. Both 
effects lead to an accelerated time dependence of the motion. Including the inertial term in the Langevin equation 
Eq. (^) changes the solution for short times according to 

m7 \ 7 

+ j dsri(s) (l - e - 7( *- s) ) . (19) 

At times much greater than the characteristic time 7 , we have 

F 



x(t) ~ —t + / dsri(s) , (20) 
«7 Jo 

which follows from the behavior of the Laplace transform of the convolution integral in the second summand of Eq. 
(|l9|): ^1 — ~^Py^) ~ m * ne um it u ^ 7- Thus the effect of the inertial term is negligible for times t much 
greater than 7 . 

If the first moment exists, that is for 1 < fi < 2, we find from Eq. fll?] ) 

<*(*)> = ^. (2D 

Only for /i = 2 we have a finite second moment, and the mean square displacement becomes ([(x(t) — (x(t))] 2 ) = 2Dt, 
in agreement with Eqs. (|l^) and ([l8|). 

For the standard FPE as well as for the subdiffusive FFPE introduced in (^8|, one finds the generalized Einstein 
relation |28|j47| , ^ | (x(t))p = Fo(x 2 (t))o/(2kBT), relating the first moment in presence of the force -Fo to the second 
moment in absence of the force. For the Levy flight model defined in Eq. ([ll]), only in the Brownian limit /1 = 2 
this relation is satisfied, provided we choose the proper amplitude of the noise, i.e. take it to be thermal noise: 
D = ksT /[^/m]. Generally since we have a diverging mean square displacement, the generalized Einstein relation does 
not hold, and we have a violation of the classical fluctuation-dissipation theorem. 

IV. LINEAR FORCE AND NON-GIBBSIAN STATIONARY SOLUTION 

In the case of ordinary Brownian motion the diffusing particles can be trapped in a harmonic potential and thus 
attain an equilibrium distribution with a finite variance | 48| , |36|| . More precisely this equilibrium distribution is the 
Gibbs or Boltzmann distribution, also obtainable from maximizing the Gibbs entropy under the constraints of norm 
and energy conservation. Also for subdiffusive transport in a harmonic potential, this property is fulfilled |28| . For 
Levy flights we shall see that a stationary solution does exist, however it possesses no finite variance. This deviation 
from the Gibbs-Boltzmann equilibrium implies that Levy flights do not describe systems close to thermal equilibrium. 

For the Hookean force F(x) = — Xx, corresponding to the harmonic potential V(x) = \\x 2 , the FFPE Eq. (|) 
becomes 
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d 



W(x,t) 



xW(x,t) + DV^W{x,t). 



dt ' dx \7to 
In Fourier-space, the conjugate equation reads 

%-W(k,t) = -—k^-W{k,t) - D\k\»W(k,t), 
at 7m ok 

which can be easily solved by making a transformation of variables (applying the method of characteristics) 

-fmD\k\v r 



W[k,t) = exp 



//A 



1 - e 



-fiXtf-ym 



(22) 



(23) 



(24) 



HOE (\ _ e -Mt/7m) 5 anc i the exact solution in real 



This is still a Levy distribution, only with a different "width" D —-j 
space can again be obtained from Eq. (|^) by inserting the time-dependent width. For /i = 2 we recover the Brownian 
results, but in the general case /i < 2 a different situation arises. We always reach a stationary distribution 



Wst (ft) = exp 



7mD|fc| A 



(25) 



but with a diverging mean square. The exact stationary solution in cc-space can be given in terms of Fox's Ti-functions: 

\x\»fj,\ 



(1.1), (14) 
(l,M),(l,f) 



(26) 



leading to the asymptotic power-law behavior W st (x) ~ I>ym/(/iA|2;| 1+Al ). 

A numerical result for a simulation of a Levy flight in an harmonic potential is shown in Fig. |^. The slope is in 
good agreement with the theoretical prediction. 

We pause to mention that we could have derived the solution Eq. ( p4| ) also by means of a separation ansatz, i.e. 
assuming a particular solution of the form W n (x,t) = T n (t)<p n (x), as was discussed in Refs. f^Jio]] . Thus we arrive 
at the ordinary differential equations 



-T(t) = -X n T(t) 



\ n fn(x) + — 

ax 



X<p n (z) 

7m 



•DW„(a;) =0 



(27a) 
(27b) 



with the eigenvalue A„. The complete solution is then given by the sum W(x, t) = W n (x, t). For the time behavior 
we find the usual exponentially decaying modes, T n (t) = e~ A ™*, and for the spatial eigcnfunction we have 

Vn {k)=c n \k\ x ^ m ' x e-^^. (28) 
The eigenvalues are given by A„ = —/in, and the complete solution in wave number space is given by 



W(k,t) = e-Dr/mW/M 

urn** 



^n^—finXt/['ym,] 



n=0 



The sum converges to Eq. (24), and can be transformed back to real space by use of the Fox functions: 



m] g — finXt I [7m] 



xH. 



n=0 

1,1 
2.2 



/zA 



1 _ 2s. 1) 



(Lm),(i,§) 



(29) 



(30) 



In comparison to Risken's result for the Fokker-Planck equation in a harmonic potential in the Brownian case [ p6[ , 
also referred to as the Ornstein-Uhlenbeck process, the eigenvalues in the solution Eq. (|30|) for \i = 2 take on only 
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even numbers. This is due to our consideration of the start in the origin, so that all the uneven Hermite polynomials 
occurring in the solution given in Ref. |36f] , vanish: i?2n+i(Q) = 0. We note that the Fox functions in Eq. ( |30|) can 
be considered as generalized Hermite polynomials. Clearly, the stationary solution corresponding to Ao = obtained 
from Eq. (30) is the same stable law as in Eq. (^6|). Also it can be seen, that only in the Brownian case [i = 2 do we 
recover the Boltzmann distribution W(x) oc e~ Xx /[ 2fc s T l. Thus Boltzmann equilibrium with a finite variance is not 
reached, in spite of the fact that the system is isolated and time-independent in respect to the ensemble. This can 
be physically understood as a consequence of the diverging mean kinetic energy of the free Levy flight. In this case 
we have 



— v (t) = -jv(t) + 777(f). 



(31) 



If we take the noise to be white according to Eq. (||) with 71 = 2, we get from Eq. (31) that in the stationary state 

rwfD 



(v 2 ) = Dj ^ (E km ) = 



(32) 



When the external harmonic potential is turned on, a length scale is introduced by the comparison (Ekin) ~ (E po t) = 
^X{x 2 ). This means that {x 2 ) « mD"//X. In fact, solving Eq. (|J) in the Brownian case with the harmonic potential, 
we obtain exactly (x 2 ) = mD^/X, in accordance with the equipartition theorem. Similar considerations remain 
valid for the subdiffusive FFPE from Ref. 45 1 , which thus describes anomalous systems close to thermal equilibrium. 
However, in the Levy case we have 



(33) 



and therefore the length scale which appears when the potential is introduced, is also diverging. Consequently the 
question arises whether other statistics could predict the equilibrium distribution in the present context. We here 
consider the recently proposed Tsallis's g-statistics [p9|, according to which the generalized entropy 



S g \p(x,v)] = 



is introduced along with the generalized constraints 



1 — J p q (x, v) dxdv 



p(x, v) dx = 1 
p q (x, v)E(x, v) dx = U: 



(34) 

(35a) 
(35b) 



for q — ► 1, S a re covers the usual Boltzmann entropy. Here v is the velocity of the particle and E(x,v) is its energy. 
Thus Eq. (35b) is a generalized constraint of conservation of energy along with the usual norm conservation Eq. 
( |35a| ). Varying Eq. ( |34| ) subject to these constraints by introducing Lagrange-multipliers, one obtains the stationary 
distribution 



p q (x,v) ~ ( 1 - (1 - q)^ [Xx 



2 1 mv 2 ] 



1/(1-9) 



Integrating over all velocities to obtain the distribution of positions alone yields 



l>,(-r) - ( 1 - (1 — q)^Xx 2 



3-9 

2-2q 



(36) 



(37) 



compare to Ref. Matching this expression to the asymptotic behavior of the stationary solution of the FFPE, 

Eq. (|25|), we obtain 



M = 



4-2g 
9-1 



(38) 



The relationship Eq. ( pq ) implies that q can range in the interval (1,2). This is at variance with the relation 
(i = (3 — <7free)/((?frcc — ljfound in the case of the free Levy flight 30 3j|, where the allowed range is (7f rce € (5/3, 3). 



G 



Note that the q-\x relation does not involve the potential strength A. Thus, the Tsallis index q makes a jump for a 
given Levy index /x, when the potential is switched off, irrespectively of how slowly, e.g. quasi-statically, the limit 
A — > is performed. Moreover, away from the asymptotic regime, the solution predicted by the entropy Eq. (|34|) does 
not agree with the solution found in Eq. (|25|) in the stationary state. Thus we conclude, that the Tsallis entropy is 
not the appropriate framework for Levy flights in a harmonic potential described by the generalized Fokker Planck 
equation Eq. ([|) . This form of a generalized entropy does not give rise to the solution of Eq. . Recently, it has 
also been shown that Tsallis g-statistics is by no means unique ]49[] , so that it would have been rather surprising if 
this very special form of statistics had led to the complete description according to the Levy flight model. Finally, 
the FFPE Eq. (^) being linear, is not compatible with the non-extensive nature of Tsallis entropy p9| ; compare with 
the non-linear diffusion equation derived from Tsallis entropy in Ref. ]5C|| . 

By comparing the distribution of the particle in the harmonic potential W(x, t) with that of the free flight Wq(x, t), 
we obtain the correspondence 

W(x,t)=W (x,t eS ), (39) 

where we define 



^ (l - e -* A */h™]). (40) 



tcff = , 

/iA 

Thus the distribution of the particle position in an harmonic potential can be obtained from the distribution in the 
free Levy flight case at an earlier, "effective" time i e fj. This comparison illustrates the slowing down of the particle 
in the harmonic potential, where the restoring force is centered towards the origin. It characterizes in a precise way 
the approach to stationarity, which is graphed in Fig. |(| Furthermore, this approach is seen to take longer time the 
smaller [i is; the quickest relaxation occurs in the Brownian case fj, = 2. 

We conclude this section by some remarks about the inertial term encountered for the Hookean force. There are 
two time scales (decay times) involved in this case: 

/ a \ \ 

(41) 




Considering only the case of large over damping, i.e. j 2 ^> 4A/m, the two time scales separate into a fast and a slow 
decaying mode, 77 — > 7 -1 and r s — > jm/X, with 77 <C t s . It can be easily shown, that neglecting the fast mode for 
times long compared to 77, corresponds to neglecting the inertial term in the original Langevin equation. Thus we 
have an effective separation into three different regimes, and being in the second with 77 <C t < r s , the approach to 
the stationary state is not influenced by the omission of the inertial term. 



V. SOLUTION OF THE LANGEVIN EQUATION 

All of the above results could have been reached equally well directly from the Langevin equation Eq. (Q). To 
illustrate this, we solve this equation for a Levy flight in a constant force field in addition to a linear force, F(x) = 
— Xx + Fq corresponding for instance to a harmonic potential and a superimposed gravity field. However, it could also 
correspond to many harmonic oscillators placed at different positions, i.e. 

N I N \ N 

F(x) = ~^i( x - x i) = - I 51 ^ J X + X! ^ lXi 
1=1 \i=l / 1=1 

= -\x + Fq. (43) 

The solution of Eq. (||) is given by 

x(t) = e^'/h™] f dt'e xt '^ ( M) + —) . (44) 
Jo V 7«V 

The distribution can then be found using the identity 



7 



P (x,t) = (S(x-x(t))) 
dk 

(exp(ik[x — x(t)])) 



(27T) 

dk 

W) 



e lkx p{k,t). 



Using the solution Eq. ( f44| ) to obtain for the characteristic function p(k, i) 

p(k,t) = {exp(-ie~ xt /^k 



x / dt'e xt '^ m] 
Jo 



nit') + ^ 
7m 



we have after discretizing the integral 

p{k,t) ~ e 



-ifeF /A(l-e- At /h m l) 



(Yl cxp f-te-^-^/h^AJtyCO)) 



(45) 



(46a) 



(46b) 



t'=0 



p(fc,t)^e-^/ A ( 1 -^ At/hH ) 

x Y[ (cxp (-ze- A(t ~*' )/[7ml Afcr7(t'))} 



Using the definition of Levy noise fromEq. (|), we obtain 



p(k,t) ~ e 



-ifeF /A(l- e - A *''[T m )) 



exp UDe-^ t - t "> /bm] A f "\k\^ 



t'=0 



(46c) 



(46d) 



Reintroducing the integrals and using the same renormalization DA^ 1 — > Z? as in passing from Eq. (||) to Eq. (Q), 
we finally have 



p(M) = e-^/^ 1 "^™) 



x exp —D'ym 



1 - e 



- fj,\t f[ym] 



\k\f 



(47) 



For A = 0, we recover the constant force result, Eq. (|l7j). For F — we get Eq. (p4|). The free Levy flight result 
Wo(x,t) according to Eq. (||) is likewise reproduced, when A = Fq = 0. In fact, by comparing to Eq. (|J) we have the 
correspondence 



W KFo {x,t) 

F r 



Wo x 



A 



1-e 



-At/b™] 



7m 
/iA 



1 - e 



(48) 



In the presence of the harmonic potential we can no longer simply make a Galilean transformation in order to eliminate 
the constant force, since the presence of the linear force singles out another special reference frame. 



VI. SOME REMARKS ON THE NUMERICAL SIMULATIONS 



Using a computer code written in C, we have simulated Levy flights in two dimensions in order to compare with 
the theoretical predictions. The noise has been defined by the asymptotics of the Levy distribution to 
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(49) 



We have basically investigated three properties: 1) histograms of the distribution of position for the free walker, 2) 
histograms of the walker in a harmonic potential, and 3) the dynamic exponent as defined via the imaginary box in 
Eq. @. 

The imaginary box according to Eq. fXLti) grows in time like the characteristic width of the stable distribution, 
{x 2 {t))i J ~ i 2 /^, which is not the variance. It gives a measure, that a finite portion of the probability is gathered 
within a given interval, which we call the imaginary box. The values of L\ and Li have been chosen so as to 1) ensure 
that we are in the asymptotic regime, where W(x, t) ~ and 2) to produce good statistics for all values of /i. 

When fitting the results to a straight line on a log-log plot, we have selected a subset of equidistant points from the 
entire set of data, in order not to favor the high t region over the low t region. 

Concerning the histogram, several precautions have to be taken when working with power law statistics. First of 
all, due to the occurrence of arbitrarily long steps, we have to define the interval of sampling beforehand, since this is 
the only way to improve the statistics when increasing the number of samples. We have chosen a minimum limit for 
the evaluated data points to ensure the asymptotic range. The maximum limit has been chosen as the maximum of 
the first say 100 (out of a total of 10000) simulations. In this way we obtain many data points throughout the entire 
region, where the asymptotic power-law expression is valid. 

We close with a remark on the axes of the plots of the histograms. To obtain equidistant points on the log-log 
plots, we chose to graph the histograms of the distributions p{y) of y, where y = log a;. If now p(x) ~ we 
have p(y) ~ e~^ y . Plotting the logarithm of p(y) as a function of y, we find a straight line with slope — fi, which is 
not to be confused with the slope — 1 — p, for the log-log plot of p(x). 

VII. CONCLUSIONS 

We have investigated Levy flights under the influence of external force fields. Especially for the cases of free flights, 
a constant and a linear force explicit solutions are derived and the consequences shown. The solutions can be reduced 
to a transformation of variables in the free flight result. We have employed the approaches of FFPEs and a Langevin 
equation with a power law noise term. We have shown that the classical fluctuation-dissipation theorem is violated in 
the case of constant force, and exhibited the non-Gibbsian nature of the "equilibrium" distribution in the harmonic 
potential. The approach to stationarity has been characterized by an effective time, and the connection to Tsallis's 
g-statistics has been explored in some detail. For the constant force, due to the divergence of the mean square 
displacement, the generalized Einstein relation breaks down, so that the FFPE Eq. (|j) also violates linear response. 
Numerical simulations have been found to agree well with the theoretical analysis. 

The direct solution of the Langevin equation Eq. (||) was given, in agreement with the solutions for the FFPE 
approach. The Langevin approach contains a priori more information than the corresponding FFPE, and might 
be more suitable for a physical interpretation of the underlying system. For a given problem, however, it will be 
more convenient to employ the FFPE approach and use the methods of characteristics or the separation of variables. 
Especially, the extraction of moments (x n ) is straightforward using the FFPE, by noting that (x n ) = J dxW(x, t)x n . 

Levy flights are typical for systems off thermal equilibrium, where thermal equilibrium is to be understood in the 
classical Boltzmann-Gibbs sense. For the systems analyzed in Refs. fl6|-|l9f|, or similar situations, the remoteness 
from thermal equilibrium is evident. For the bulk mediated surface diffusion process in Ref. one should keep 

in mind, that the Levy process emanates only as an effective motion, resulting from a Brownian walker which is 
eventually absorbed on the surface, before he gets activated again. 

The connections to the Tsallis's g-entropy still remain unclear. The discrepancy of the relations of q to the Levy 
index /i in the free and harmonic cases show, that the concept of g-entropy cannot provide a full explanation of Levy 
flights. Especially the stationary solution found for the harmonic potential, is not an equilibrium solution in the sense 
of Tsallis's q-entropy. To our understanding, only a non-linear generalized FPE can lead to results compatible with 
this theory. 

We believe that our analysis provides further understanding of anomalous diffusion processes, and will give rise to 
further experimental investigations, for example Levy type reaction dynamics subject to an electric field, or the tracer 
diffusion in rock structures under gravitation, or similar. 
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APPENDIX A: CORRELATION FUNCTIONS 



Here we examine two-point correlations like (x(t)x(t')) in the harmonic case, according to the tools developed in 
Sec. [v|. They are often sufficient for the description of the system in applications of the underlying (F)FPE. Assuming 
the initial condition a;(0) = 0, we introduce the Green's function 

(Al) 



(A2) 



x(t) = J dt' G(t,t')r](t'), 
and by comparison with Eq. (^4j) it is seen that 

G(t,t') = e^-^/^e^ - t')e(i'). 



Due to the divergence of the moments of the Levy distribution in general, we have to work with characteristic functions 
discussing correlations: 

/g-» J dsA(s)x(s)^ 

= (exp(-i [ ds [ ds' A(s)G(s,s')t){s'))) 



(A3) 



TT^ e -iAj7( S ') J dsA(s)G(s,s')^ 
TT--D|A J dsA{s)G(s,s')^ 



where A(t) denotes an a priori arbitrary function, and we used the definition Eq. (g). With the usual renormalization 
of the noise DA M_1 -> D |E1 we have 



^g-» / dsA{s)x(s)^ 



exp (-D J ds' J 
exp J ds' 



dsA(s)G(s,s') 



ds J ds"G(s,s')G(s" 7 s')A(s)A{s") 
exp ( -D / e^'^W 



ds 



x / ds"e- x{s+s " )/[ ~< m] A{s)A(s") 



In the last step we used the identity 

G(S,S')G(S", S ') = e -Ms+s"-2s')/[ r m] 

xQ(s-s')Q(s"-s')Q(s'). 

Putting A(s) = A8(t - s), we find 



(A4) 



(A5) 
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-iAx(t)\ 



exp —D 



exp 



x e 



exp [ — D 
exp I DA 1 ' 



ds ' e H\s'/[ im ] J dg 

da" e -^'+'"y^A 2 5(a - t)6(s" 
ds' 

e -2\t/hm] A 2 e( j._ s > 
^ se -M(t-s)/[7m]^At 

7 m(l - e -Mt/b™-]) 



/a As /[7m] 



M/2 



//A 



(A6) 



which is equivalent to Eq. (]24|), as it should be, and therefore includes also the results for the Brownian case. One- 
point correlation functions (moments of the distribution), if they exist, can thus be obtained from Eq. (A6), so we will 
proceed to the more in teres ting case of the two-point correlations. To this end, we take A(t) = A [5(t — tij — S(t — t 2 )] 
and insert it into Eq. (A4), and we find 

(exp (-aO(ii) -x(t 2 )))) = 

" 7 m(l - e-^/M) 



exp 



-DA» 



1 - e 
7m(l — e 



//,A 

fj,X(ti— t2)/[jm] 



//A 



(A7) 



for t\ > and t\ and i 2 interchanged in the other case when i 2 > t±. This is essentially the characteristic function of 
the stochastic va riab le x{t\) — a;(i 2 ), so all two-point correlation functions and the distribution itself can in principle be 
found from Eq. (A7). However, it is seen that it is a Levy distribution, so all higher moments diverge. Nevertheless, 
Eq. (A7) still gives some information about the correlation between the position of the walker at two different times. 
When ti i 2 the characteristic function splits up into the product of the characteristic function of the two variables 
Eq. (A6) which means, that x(t\) and a;(i 2 ) are independent in this limit. At intermediate times, i.e. when both 
t\ and i 2 are small, the correlation depends on both. This is a memory of the initial conditions, since both walkers 
start out at the origin. At long times, initial conditions are not important, and hence only a dependence on the time 
difference t\ — i 2 (kept fixed and finite) is retained, 



-t^(x(ti)-x(t 2 )) 



exp 



) 

-Djm 
/iA 



-fi\\t 1 -t 2 \/[ym] 



_l _ g-MI*i-*2|/[7ra] 



Writing X12 — x(ti) — a;(t 2 ), we have by the usual arguments a Levy distribution of cci 2 : 



(A8) 



(A9) 



with 



7m 



D M (*i,t 2 ) = D^- 
/iA 



-li\t 2 /[-ym] ^ _ e -p.\(t 1 -t 2 )/[-ym]y 



-/iA(ii-t 2 )/[7m] 



(A10) 



As with the other Levy distributions, this stochastic variable is characterized by a power law tail (ti = t, t% = 0): 
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p(x 12 ,t) ~ D-^- ((l-e-^l^y 

+1 _ e ^t/hm]"j | Xl2 |-l-f. (All) 

In the case of /i = 2 all the preceding results reproduce the well known Brownian relations for ordinary diffusion. 
One could proceed like this finding three-point correlations in an analogous manner. 
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FIG. 1. Typical Levy-flight for the Levy index [i — 1.4. The clustering is obvious. Each cluster is statistically self-similar 
to the unmagnified picture. The fractal dimension of the flight is df = /J, plM- 
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FIG. 2. Function (x 2 (t))L from Eq. ( [il"|) versus time with fj, = 1 in a log-log plot. The slope of the straight line is 
1.990 ± 0.028, which is to be compared to the expected value 2//i = 2. 
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FIG. 3. log-log plot of {x 2 ) L versus time from Eq. (jH]) (lower curve), compared to the squared absolute mean (\x\) 2 from 
Eq. ( p^ ) (upper curve) for n — 1.5. The fitted slope is 1.363 ±0.013 and 1.364 ±0.035, respectively, which is in good agreement 
with the theoretical value 4/3. 
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FIG. 4. Graph of the slope 2/fi according to Eq. ( |ll| ) as a function of the Levy index [i. Note the bend at /i = 2 marking 
the transition to normal diffusion. 
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FIG. 5. Histogram for the stationary solution W s t of a Levy flight in a harmonic potential versus |as | , as a plot of log W s t(y) 
versusj/, where y — log(x) denotes the natural logarithm of the position of the flyer, see Sec. VI. The data were produced for a 



Levy index \i = 1.4. The fit indicated by the dashed line reveals a slope of —1.408 ± 0.108, which thus shows a good agreement 
with the theoretical prediction — fj,. 
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FIG. 6. The linear time as seen by the free walker (dashed line) , compared to the effective time sensed by the random walker 
in the harmonic potential, as a function of the laboratory time. For the Levy flight in the potential, the restoring Hookean 
force slows down the spreading of the diffusion, and eventually brings it to a halt. 
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